wd<- "/Users/chengg/Google Drive/EPICON/Mycobiome/Fungal ITS/statistic/Total.fungi"
setwd(wd)
rm(list = ls())
load("EPICON.data.preparation.RC.bNTI.ted.2019.04.19.Rdata")
library(vegan)
library(ggplot2)
library(colorRamps)
library(ape)
library(umap)
library(plyr)
library(plotly)
set.seed(315)
fung.pc<-pcoa(vegdist(decostand(fung.rar,"hellinger"),"bray"))
envpc<-data.frame(fung.pc$vectors, env)
pc1=round(fung.pc$values$Relative_eig[1], 3)*100
pc2=round(fung.pc$values$Relative_eig[2], 3)*100
envpc$Habitat<-factor(envpc$Habitat, levels = c("Leaf", "Root", "Rhizosphere", "Soil") )
envpc$Treatment1<-factor(envpc$Treatment1, levels = c("Control", "Pre_flowering_drought", "Post_flowering_drought"), labels = c("CON", "PRE", "POST") )
df <- envpc[,c("Axis.1", "Axis.2", "Habitat")]
find_hull <- function(envpc) envpc[chull(envpc$Axis.1, envpc$Axis.2), ]
hulls <- ddply(df, "Habitat", find_hull)
p1<-ggplot(data=envpc, aes(x= Axis.1 , y= Axis.2,group=Habitat))+
geom_point(aes(shape=Habitat, colour=factor(TP)),size=2) +
geom_polygon(data=hulls, alpha=.1)+
labs(x = sprintf("PCo1 (%.1f%%)", pc1), y = sprintf("PCo2 (%.1f%%)", pc2))+
scale_shape_manual(values=c(0, 6, 3, 1))+
scale_colour_manual(values=blue2red(18))+theme_bw()+
guides(color=guide_legend(title= "Week"),
shape=guide_legend(title= "Compartment"))+
theme(legend.title = element_text(colour="black", size=15, face="bold"),
legend.text = element_text(colour="black", size=15, face="bold"),
axis.text=element_text(size=10,face="bold"),
axis.title=element_text(size=15,face="bold"))
p1

ggplotly(p1)
library(plotly)
plot_ly(x=envpc$Axis.1, y=envpc$Axis.2, z=envpc$Axis.3, type="scatter3d", mode="markers", color=envpc$Habitat)
envpc.sub<-envpc[envpc$TP<2 &envpc$Treatment1=="CON",]
df <- envpc.sub[,c("Axis.1", "Axis.2", "Timepiont", "Habitat")]
df$TPHabitat<-droplevels(interaction(df$Timepiont, df$Habitat))
find_hull <- function(envpc.sub) envpc.sub[chull(envpc.sub$Axis.1, envpc.sub$Axis.2), ]
hulls <- ddply(df, "TPHabitat", find_hull)
p2<-ggplot(data=df, aes(x= Axis.1 , y= Axis.2,group=Habitat))+
geom_point(aes(shape=Habitat, colour=Timepiont),size=2) +
geom_polygon(data=hulls, alpha=.1)+
labs(x = "PCo1", y = "PCo2")+
scale_shape_manual(values=c(3,17,16, 15))+
scale_colour_manual(values=c("blue", "red", "brown"))+theme_bw()+
guides(color=guide_legend(title= "Week"),
shape=guide_legend(title= "Compartment"))+
theme(legend.title = element_text(colour="black", size=15, face="bold"),
legend.text = element_text(colour="black", size=15, face="bold"),
axis.text=element_text(size=10,face="bold"),
axis.title=element_text(size=15,face="bold"))
p2

ggplotly(p2)
subs<-((envpc$TP == 8 & envpc$Treatment1 != "POST") |
(envpc$TP == 17 & envpc$Treatment1 != "PRE")) &
(envpc$Habitat=="Root" | envpc$Habitat =="Rhizosphere")
p3<-ggplot() + geom_point(data=envpc[subs,],aes(x= Axis.1 , y= Axis.2,shape=factor(TP), colour=Treatment1),size=3,alpha = 0.6) +
labs(x = sprintf("PCo1"), y = sprintf("PCo2"))+
scale_colour_manual(values=c("black", "blue", "red"))+theme_bw()+
facet_wrap(~Habitat,scales="free", nrow=1)+
guides(color=guide_legend(title= "Treatment"),
shape=guide_legend(title= "Week"))+
theme(strip.text.x = element_text(size = 15,face="bold"),
legend.title = element_text(colour="black", size=15, face="bold"),
legend.text = element_text(colour="black", size=15, face="bold"),
axis.text=element_text(size=10,face="bold"),
axis.title=element_text(size=15,face="bold"))
ggplotly(p3)
set.seed(315)
fung.umap<-umap(decostand(fung.rar,"hellinger"))
envumap<-data.frame(fung.umap$layout, env)
envumap$Habitat<-factor(envumap$Habitat, levels = c("Leaf", "Root", "Rhizosphere", "Soil") )
envumap$Treatment1<-factor(envumap$Treatment1, levels = c("Control", "Pre_flowering_drought", "Post_flowering_drought"), labels = c("CON", "PRE", "POST") )
p4<-ggplot() +
geom_point(data=envumap,aes(x= X1, y= X2, shape=Habitat, colour=factor(TP)),size=2)+
labs(x = "Umap1", y = "Umap2")+
scale_shape_manual(values=c(0, 6, 3, 1))+
scale_colour_manual(values=blue2red(18))+theme_bw()+
guides(color=guide_legend(title= "Week"),
shape=guide_legend(title= "Compartment"))+
theme(legend.title = element_text(colour="black", size=15, face="bold"),
legend.text = element_text(colour="black", size=15, face="bold"),
axis.text=element_text(size=10,face="bold"),
axis.title=element_text(size=15,face="bold"))
p4

ggplotly(p4)